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Abstract 

We determine the strong coupling constant a s from a lattice calculation of vacuum polarization 
functions (VPF) in three-flavor QCD with dynamical overlap fermions. Fitting lattice data of VPF 
to the continuum perturbative formula including the operator product expansion, we extract the 
QCD scale parameter A^L. At the Z boson mass scale, we obtain a^(Mz) = 0.1181(3)(ij2)> 
where the first error is statistical and the second is our estimate of various systematic uncertainties. 
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I. INTRODUCTION 



Strong coupling constant a s is one of the fundamental parameters in the Standard Model. 
Its precise determination from various sources provides a crucial test of Quantum Chromo- 
dynamics (QCD). Experimentally, it is obtained through high energy particle scatterings 
involving quarks, such as e + e~ — > hadrons (see [l, 2] for a summary), for which perturbative 
calculation of QCD is possible. 

Among other measurements, the hadronic decay rate of tau lepton jsj provides one of 
the most precise determinations of a s . The tau- lepton hadronic decay rate is written in 
terms of a vacuum polarization function (VPF) of weak currents. Since the perturbative 
QCD calculation cannot be directly applied for physical time-like momentum transfer (or 
the virtual W boson mass) q 2 , one considers a total decay rate that involves an integral 
over q 2 . The decay rate is then calculated using three- loop perturbative expansion and the 
Operator Product Expansion (OPE) to parametrize corrections at low energies. The final 
error in a s contains those from the perturbative expansion and uncertainties of condensates 
appearing in OPE, both of which are non-negligible because the relevant energy scale is of 
order of 1 GeV, i.e. lower than the tau lepton mass m T . There is also an assumption of the 
quark-hadron duality, which is not trivially satisfied. 

Theoretically, the perturbative calculation and the extraction of a s become more trans- 
parent if the experimental data for VPF are available at space-like momenta. Although there 
is no such direct measurement, lattice QCD is able to provide a non-perturbative calculation 
of VPF at space-like momenta. Since the calculation is based on the first-principles of QCD, 
for which the input parameters are low-lying hadron spectrum, it gives another method to 
extract a s . 

Using lattice QCD, one can calculate VPF of vector (V) and axial-vector (A) channels 
from two-point functions of those currents. After a Fourier transformation in four space-time 
dimensions, VPF can be obtained as a function of space-like momentum squared Q 2 ranging 
from zero to the order of lattice cutoff squared. Our recent lattice study with two flavors 
of dynamical overlap fermions demonstrated that such data can indeed be used to extract 
a s combining with a perturbative calculation in the continuum scheme ^|. The result is 
consistent with other lattice calculations in two-flavor QCD {5, 6]. 

One of the main advantages of our method is that the calculation can be done on existing 
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gauge configurations produced for light-hadron spectrum Furthermore, unlike other 
methods previously applied, there is no need of multi-loop perturbative calculation on the 
lattice, which is so complicated that one typically has to develop automated tool dedicated 
for a given lattice action. Finite volume effect is negligible for VPF at relatively large Q 2 , 
and discretization effect is carefully studied and is shown to be under control. The method 
therefore provides a reliable and economical way to extract the most important parameter of 
QCD, i.e. a s . A similar idea has also been applied for the charmonium two-point function 
to determine a s as well as charm quark mass [9|. 

In this work, we extend our previous study to the case of realistic three-flavor QCD 
with dynamical light and strange quarks. We also improve the calculation by employing 



the conserved current for the overlap fermion 



lOj , which simplifies the possible form of the 



two-point function as it satisfies the Ward-Takahashi (WT) identity on the lattice. The 
extraction of VPF is thus made more straightforward. The determination of the strong 
coupling constant uses the continuum perturbative QCD formula up to four loops and OPE 
up to 1/Q terms. The value of the quark condensate is calculated independently and used 
as an input in this work. The result is translated to the common definition, i.e. a^(Mz), 
the strong coupling constant of five-flavor QCD in the MS scheme at the Z boson mass scale. 
Numerical simulations of lattice QCD are carried out with 2+1-flavors of dynamical 



fermions described by the overlap fermion formulation [111, [12]. We have data at a single 
lattice spacing with the inverse lattice spacing a -1 = 1.83(1) GeV estimated using the static 
quark potential with an input for the Sommer scale r = 0.49 fm (the associated error for this 
is discussed later). The lattice volume is 16 3 x 48, which leads to the physical volume about 
(1.8 fm) 3 x(5.4 fm). The gauge configurations are generated in the course of the dynamical 
overlap fermion simulations by the JLQCD and TWQCD collaborations The set of up 
and down quark masses m u d covers the range of (0.2-0.8)mP hys and the set of the strange 
quark mass m s covers the range of (1.0-1.3)mP hys with mf ys the physical strange quark 
mass. The valence quark mass in the calculation of VPF is set equal to the up and down sea 
quark mass. For each combination of m u d and m s , we use 260 configurations, each of which 
is separated by 100 HMC trajectories. Global topological charge of the gauge configurations 
is fixed to zero, which may induce small finite size effect for long distance physical quantities 
[ijj]. We expect that this gives negligible effects on the short distance physics considered in 
this work. 
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The rest of the paper is organized as follows. In Section |TT1 we describe the details of the 
lattice calculation of VPF, including the definition of the overlap fermion used in this work 
and a construction of the conserved current for the overlap fermion. Section II III discusses 
the fit of VPF to the perturbative formula. The possible systematic errors are discussed in 
Section IIV} followed by our final results. Our conclusions are given in Section |Vj 



II. LATTICE CALCULATION OF VACUUM POLARIZATION FUNCTIONS 

In the continuum theory, transverse (T^j\Q)) and longitudinal (IIj (Q)) parts of VPF are 
defined through two-point functions (J^(x)J^(0)) of either vector (J^ = V^) or axial-vector 
(J M = A^) currents with a, b the flavor indices. Namely, after a Fourier transformation to 
the momentum space, the two-point functions are parametrized as 

(j;4)(Q) = S ab [(5, U Q 2 - Q,Q u )Ilf{Q) - Q^ufiQ)] , (1) 

where the momentum is space-like as we work on an Euclidean space-time lattice. Be- 
cause of the WT identities, the longitudinal part of the vector channel vanishes, Hy\Q) = 0, 
while the axial-vector channel is proportional to the quark mass. 

On the lattice, we employ the overlap fermion formulation [ljj], whose action S ov = 
'}2 xy q(x)D ov (x,y)q(y) is specified by the massive overlap-Dirac operator 



D ov (x, y) = (m + —J + (m - yj l5Sgn[H w (x, y; -m )] (2) 

for a quark mass m. Here, m is a parameter to define the overlap kernel Hyy(x, y; —m ) = 
, y 5 (Dw( K x,y) — mo) with Dw(x, y) the conventional Wilson-Dirac operator. In this study 
we take m = 1.6. (Here and in the following we set a = 1, unless otherwise stated.) The 
overlap action has an exact symmetry under a chiral rotation defined with the modified 
chirality operator 7 5 (x,y) = ^5(S Xty — D ov (x, y)/mo), so that the continuum-like axial WT 
identities are hold on the lattice at finite lattice spacings. 

The conserved vector current for this action has a complicated form, which can be written 
in a general form V® ,cv (x) = z q{w)K lJi {w, z\x)q{z) with a non-local kernel K^w, z\x). 
K^w, z\x) is determined such that it forms a Noether current under a local vector transfor- 
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mation 

5S ov = ^2q(x)[-a(x)D ov (x,y) + D ov (x,y)a(y)]q(y) 

x,y 

= J2^M^K,(x,y\z)q(y) (3) 

x,y,z 

with a(x) a local parameter The derivative d** denotes a backward derivative operator 
d^*Vfj,(x) = V^(x) — V^(x — aft) on the lattice. Similarly, for flavor non-singlet transforma- 
tions, we can derive flavor non-singlet conserved vector and axial-vector currents as 

V*' cv (x) = Y,^yK,{w,z\x)q{z), (4) 

w,z 

A a ^(x) = Y,QHt a K,(w,z\x)[%q}(z), (5) 

w,z 

where t a denotes the generator of SU(Nf) normalized as Trt a t b = 5 ab /2. For flavor SU(2), 
t a = T a /2 with r 2 the Pauli matrix. In the following, we consider the flavor non-singlet 
currents. 

In practical implementation for numerical calculations, we approximate the sign function 
in (T5]) by a rational function with Zolotarev's optimized coefficients. In our setup, the sign 
function is approximated to the level of lCT 1 - 7 " 8 ) with the number of pole N ~ 10. (For 



details, see 



7| for instance.) Accordingly, the kernel K^(w, z\x) is constructed as 



K^(w, z\x) = m ( 1 - 7^- ) 75 



in 



a N 



-W^(h 2 w + C 2n ) 



1=1 



d N b 



i=i 

with 



A--VV / \ — — J. ~£ILJ j '2 . \" ' fJi " VV 1 - " VV ' ' fJr / J ^ I 

min 'L'W 
+ C 2 l-1 (l w + C 2 l-1 



(6) 



W^(z,w\x) = -^75 {(1 - r r fl )U IM (z)5 x+fljW 5 XtZ - (1 + i^)Ul(z - p)8 XtW 8 z ^ x ) , (7) 

and hw(w, z) = H w (w, z\ — m )/X niin , where A min is a lower limit of the eigenvalue of \H W \ 
to be approximated by the rational function. (In ([6]) the site indices of Wp and hw are 
omitted, but they are multiplied as matrices.) The Zolotarev's coefficients bi, d$ are 
given in 

In this study, we consider the two-point functions of flavor non-singlet conserved and local 
currents ( J"' cv (x) J£' loc (0)), where J^ loc {x) is either V^ loc {x) = Zq(x)t a llJL q(x) or A a > loc (x) = 



Zq{x)t ar )^^,q{x). The renormalization constant Z needed for the local currents to match 
their continuum counterpart is non-perturbatively determined as Z = 1.39360(48) flol ]. 
The lattice calculation of (J^ ,cv (x) J„' loc {y)) is standard except for the complicated form 
of J^ ,cv (x). Namely, we calculate the quark propagator originating from a fixed space-time 
point y = and construct the two-point function with the conserved current J®' cv (x) located 
at arbitrary space-time point x. We then apply the Fourier transform in all four dimensions 
to obtain the two-point function in the momentum space. 

Because of the current conservation of J^ ,cv , we may derive the WT identities for the 
two-point functions 

E^(C C X Moc )(Q)=0, (8) 
Y,QMT Ab u° C ){Q) - 2m q (P a A b / oc )(Q) = 0, (9) 

where aQ^ = sin(aQ At ) are a momentum definition corresponding to the backward derivative 
operator d**. We use a convention that the two-point function after the Fourier transfor- 
mation, such as (J^Jt)(Q): is a function of aQ^ = l-an^jL^ with £^=1^4 the extent of the 
lattice in the /x-th direction. The second term in (Q represents the correlation function of 
the pseudo-scalar density operator P a (x) = q(x)t a ^(l — D ov /mo)q(x) and the local axial- 
vector current A h J Loc {y). A possible term arising from the axial transformation of J„ loc (y) 
(J = V or A) vanishes when we take the vacuum expectation value, since the vacuum has 
axis- interchange symmetry while the index v remains in J„ loc {y)- 
The vector and axial-vector VPFs are now given by 

(J a, CVJ Moc )(g) = S a b j ((W Q2 _ Q {Q) _ Q^fiQ) + Aj,(Q)] . (10) 

Here, Uy\Q) vanishes because of the conservation of V^' cv , while IT^Q) represents a 
remnant due to PCAC: 

uf(Q) = -2m q (P a A a / oc )(Q)/(Q 2 Qu). (11) 
(Repeated indices a's are not summed.) The transverse part Hj (Q) can be extracted as 

n^(Q) = (J a / vja , loc )(Q)/(Q 2 ~ Qi*Qi*)> (12) 

(repeated indices /i's are not summed) if one ignores the additional term A^(Q), which 
reflects the violation of the current conservation of the local current J" ,loc . Since the current 
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conservation is recovered in the continuum limit, this term can be expanded in terms of 
small aQ^ as 

= £ (s,^Ql m -Ql (m ^Q^)Ql nF ^n(Q)- (is) 

m,n=l p 

where F mn denotes the scalar function depends on the index m, n and momentum Q. It 
satisfies the condition ^ Q^A^^Q) = coming from the WT identity for J®' cv . In this 
work, we confirmed that this term is numerically negligible in the range {aQ) 2 < 1, and 
ignore its contribution as we discuss later. 



III. FIT WITH THE PERTURB ATIVE FORMULA 

Defining Uj(Q) = uf{Q) + U { J\Q), the Operator Product Expansion (OPE) of VPF, 
U V+A {Q) = U V (Q) + Il A (Q), is given by 

rV+A|oPE(Q 2 ,a s ) = c + C (Q 2 ,n 2 ,a s ) 

+ Cl +A (Q 2 ,» 2 ,a s f 2{Q) 



Q 



2 



Q 4 

q=u,a,s 

+ c 0cW ',«.)<« + O (g-') (14) 

for large Q 2 . The perturbative expansion of the coefficients C^ +A ^ (X = 0, qq and GG) is 
known up to two- to four-loop order in the continuum renormalization scheme, i.e. the MS 
scheme, depending on the terms. 

The first term c in f TH|) is a scheme- dependent constant, divergent in the limit of infinite 
ultraviolet cutoff. For the Adler function D(Q 2 ) = —Q 2 dU(Q 2 )/dQ 2 , which is a physical 
observable, the first term disappears and the contributions from other terms remain finite. 



and third terms are perturbatively calculated to four-loop 
; the expression explicitly contains a^(Q) defined in the MS 
scheme. (The superscript (3) stands for the number of flavors.) The third term contains 



The coefficients in the seconc 
order in the MS scheme 16-18 



the running mass rh(Q) whose anomalous dimension is known to three- loop order [19|, |20 |. 
The fourth and fifth terms represent higher order effects in OPE containing dimension-four 
operators. Their Wilson coefficients are calculated at three- loop order [20I . 
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In addition to the terms represented by the continuum OPE (|T4j) . there are discretization 
effects of 0(a 2 Q 2 ) at finite lattice spacings. These can be eliminated by an extrapolation 
to the continuum limit, in principle. In our calculation obtained at a single lattice spacing, 
however, the error has to be carefully investigated. We use a lattice perturbation theory to 
estimate the discretization effects at large a 2 Q 2 regime as described below. We also note that 
the exact symmetries of the overlap fermion partly eliminate unphysical terms of 0(a 2 Q 2 ) 
that violate the WT identities [4|. We therefore use (I14p without including correction terms 
describing the discretization effects when we fit the lattice data of VPF extracted through 
(II ip and (fT2|) . In our previous study in two-flavor QCD (4], we had to use more complicated 
method to extract the physical VPFs, because of non-conserved (axial-)vector currents. 

We now discuss a fit of the lattice VPF data to the OPE formula. In this analysis the 
renormalization scale is set to \i = 2 GeV when necessary, though the final result should not 
depend on \i up to higher order perturbative corrections. The gluon condensate ((a s /n)GG) 
is defined only through the perturbative expression like ( [141) because of the renormalon 
ambiguity [21|, hence we treat ((a s /iT)GG) as a free parameter to describe associated l/Q A 
corrections. On the other hand, the quark condensate (qq) is well-defined in the massless 
limit, as there is no mixing with lower dimensional operators because of the exact chiral 
symmetry of overlap fermion. Thus, the quark mass dependence of Uv+a\ope{Q 2 , « s ), which 
comes from the third and fourth terms in (Tl4l) . is given only as a function of a s , once the 
quark condensate is determined elsewhere. 

The running quark mass rh(Q) is set to the value corresponding to the quark mass used 
in the lattice calculation. First, we obtain the value at 2 GeV using the non-perturbatively 
calculated Z-factor as m(2 GeV) = Z m (2GeV)m g with Z m (2GeV) = 0.806(12) (24) (j^) 



221 ] . Then, it is evolved to Q 2 using a three-loop running formula 19|, [20] 



For the quark condensate of up and down quarks, we use the value obtained in the recent 



analysis of the spectral density 



23|, i.e. (qq) = - [0.242(04) GeV] 3 , which is defined in 



the MS scheme at 2 GeV. The strange quark condensate (ss) appears only as a contribution 
from sea quark and the associated coefficient CY S +A (Q 2 , ct s ) starts from 0(a s ). For the value 
of {ss), we use the same value as the one of up and down quarks. 

In the fit of VPF using (fT4|) . there are three unknown parameters, a s , c and {(a s /7r)GG) . 
The QCD scale controls the running coupling constant a^(Q), which is evaluated 



using the four-loop formula [24j, |25 |. 
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FIG. 1: (aQ) 2 dependence of VPF, Hv+a(Q), at all valence quark masses: m q = 0.015 (circle), 
0.025 (square), 0.035 (diamond), and 0.050 (triangle). Top half is a result at m s = 0.08 while the 
bottom is at m s = 0.10. Solid curves show a fit function at each quark masses. Filled symbols are 
the points for which each momentum component is equal to or smaller than 27r/16 in the lattice 
unit. 
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FIG. 2: Comparison of Hv+a(Q) with different momentum definitions. Lattice data at m q = 0.015. 

Figure [I] shows a (aQ) 2 dependence of Uv+a(Q) in a window 0.4 < (aQ) 2 < 1.0. Fit 
curves shown in this plot are those of ffT4"j) with the value of parameters extracted from the 
fit in the range 0.463 < (aQ) 2 < 0.994. The upper limit of the range is chosen to avoid 
significant lattice artifact, which is estimated by a difference of the lattice momentum aQ^ 
from the other definition aQ^. In fact, the result is unchanged within la level when we use 
these different definitions of the momentum as Figure |2] shows, as far as (aQ) 2 is lower than 
1.0. Beyond this value we observe significant deviations between the different definitions. 
We also impose a constraint aQ^ < tt/4 for each momentum component to avoid large lattice 
artifacts. 

In order to determine the lower limit, we investigate the stability of the fit results. Figure [3] 
shows the dependence of fit parameters on the value of the lower limit (aQ) 2 ^^. We observe 
that around (aQ)min = 0.4-0.5 all the parameters are stable. 

It is interesting to see where the 1/Q 4 term becomes significant. In Figure [3j the fit 
results without the 1/Q 4 terms are also shown by filled symbols. It turned out that the 
value of is consistent with the 1/Q 4 fits when (aQ) 2 ^^ is greater than 0.5, which means 
that the 1/Q 4 terms become relevant below this value. In fact, if we extend this fit including 
the data points slightly below (aQ) 2 = 0.5, the value of A^L becomes significantly lower; the 
X 2 /dof of the fit becomes too large (~ 3.0), which suggests that the fit is no longer valid. 
Strictly speaking, x 2 /dof does not have a statistical meaning as the correlation among the 
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FIG. 3: Dependence of the fit parameters on the lower limit of the fit range. The maximum value 
is fixed at (aQ) 2 ~ 0.994. Open and filled symbols show the results with and without the l/Q A 
terms in (|14|) . (Thus, there is no filled symbol in the middle plot.) 



data at difference Q 2 is ignored in the fit used here. We discuss on the statistical correlations 
among the data points in the next section. 

The limitation of the OPE formula including up to l/Q 4 terms can be investigated by 
looking at its departure from the lattice data at lower values of Q 2 . In Figure H] we plot the 
difference of the lattice data from the fit curve including the l/Q 4 terms. The Q 2 region 
in this plot is extended towards the value below (aQ)^. From this plot, we observe that 
the next order l/Q 6 contribution becomes significant below (aQ) 2 ~ 0.4. We therefore set 
(a<5)min = 0.463 in our analysis including the l/Q 4 terms. 

After doing a simultaneous fit of the VPF data at all sea quark masses, the QCD parame- 
ter is obtained as A^L = 0.247(5) GeV. By matching onto four- and then to five-flavor QCD 
at charm and bottom quark masses respectively, the strong coupling constant is obtained as 
ai 5 \Mz) = 0. Il8l (3) at the Z boson mass scale. Here, the error is statistical only. Various 
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FIG. 4: Difference between the lattice data and the fit function (|14p . Dashed line shows a guiding 
line representing the 1/Q 6 behavior. 



sources of the systematic error are discussed in the next section. 

IV. SYSTEMATIC ERRORS AND FINAL RESULT 
A. Uncorrelated fit 

First of all, our fit procedure of VPF may induce systematic error due to the use of 
uncorrelated fit. Namely, in the fit described above, we did not take the correlation among 
different Q 2 points into account and estimated the statistical error for the fit parameters 
using the jackknife method. 

In order to estimate the associated error, we calculated the statistical correlation of 
different Q 2 points and found it very strong (50-100%). If we construct \ 2 taking account 
of the correlation, the value of x 2 /dof is of order 100. This unacceptably large value occurs 
because the fit function ( TT4"|) does not contain the discretization effects that violate Lorentz 
symmetry. Indeed, if we restrict the data points for those that each momentum component 
is equal to or smaller than 27r/16, the x 2 /dof is reduced to 1.7, without changing the central 
values of the fit parameters. The restricted data points are shown in Figure [1] by filled 
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symbols. 

In the main analysis we use all the data points that satisfy the condition aQ^ < 4ir/lQ in 
order to use as much information from the lattice data as possible with the uncorrelated fit. 
In particular, we can take a wider range of Q 2 with this choice, that improves the stability 
of the fit. In other words, with the limited data points (Q^ < 2tt/1Q) the \ 2 fit is sometimes 
trapped in a local minimum depending on the initial values for the fit parameters. 

We therefore decided to use the uncorrelated fit for the enlarged data points (Q^ < 47r/16) 
to obtain the fit parameters, and then to check the value of x 2 /dof for the limited data points 
{Qfi < 27r/16) taking account of the correlation. Since the value of x 2 /dof is 1.7, we do not 
expect the bias due to this procedure larger than one standard deviation, assuming that 
the full correlated fit should give x 2 /dof ~ 1. Thus, we conservatively assign a systematic 
error ±0.003 for af\Mz), which is equal to the size of the statistical error of one standard 
deviation. 

This procedure can be avoided if the lattice data are obtained at finer lattice spacings so 
that one can cover the same physical range of Q 2 with smaller lattice momenta. 

B. Discretization effect 

As described above, the discretization effect is significant in our lattice data especially 
when we try to cover large enough Q 2 range. We estimate the associated error using lattice 
perturbation theory. 

Since the discretization effect is most important in the large momentum region, the 
perturbation theory can be used to estimate its size. We calculate the one-loop diagram 
of VPF, Uy^_ A (Q), with local and conserved currents in lattice perturbation theory, and 
compare them with the continuum perturbation theory. This provides an estimate of the 
discretization effect at the zeroth order of a s . Because the discretization error itself is a 
small effect, its calculation at the leading order gives a reasonably precise estimate. 

Figure [5] shows the result as a function of (aQ) 2 . Since the lattice regularization violates 
rotational invariance, the result is not a completely smooth function of (aQ) 2 , as shown 
by squares in the plot, which correspond to representative points of (aQ) 2 in our lattice 
calculation. 

The perturbative result may be parametrized as ny +J 4(Q 2 ) latt ' pcrt = c — l/(2ir 2 ) ln(aQ) 2 + 
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FIG. 5: (aQ) 2 dependence of one-loop VPF Hj=v,a(Q 2 ) m lattice perturbation theory. Dashed line 
shows the leading logarithm term plus a constant, which corresponds to the continuum perturbation 
theory. Solid lines show the function including lattice artifact of 0((aQ 2 )). The shaded band 
represents an uncertainty due to the higher order effects. The red diamond denotes the value at 
the upper limit of our fit of VPF. 

0.0062(40) (aQ) 2 for small a. The logarithmic term is the same as in the continuum per- 
turbation theory and c is the scheme dependent constant as already noticed. The term 
+0.0062(40) (aQ) 2 comes from the discretization effect. The error includes a fluctuation of 
numerical integral as well as the non-smooth behavior due to higher order discretization 
effects. The non-smooth behavior appears because of different assignments of momentum 
components aQ^. In the plots we took several values of (aQ) 2 (and so aQ^) that also ap- 
pear in the lattice calculation. We observe that the pattern of non-smooth behavior in 
the one-loop calculation actually well reproduces that occurring in the numerical simula- 
tion. It suggests that our estimate of the discretization effect is reasonably realistic. As 
the plot shows, this error band is taken so that the result of the lattice perturbation theory 
calculation is covered. 

By subtracting this estimate of the 0((aQ) 2 ) effect from the lattice data, the final result 
for a^(Mz) changes by +0.0002(1). We therefore take +0.0003 as our estimate of the 
systematic error from this source, without changing the central value to be conservative. The 
estimated error in the negative direction is thus taken to be zero. Although the perturbative 
calculation is done only at the one-loop level, we expect that the higher order effects are 
suppressed by an additional factor of a s and thus well below ±0.0001. 
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FIG. 6: Difference between IIy ag (Q) and Tly^ d (Q) at all valence quark masses m q = 0.015 (circles), 
0.025 (squares), 0.035 (diamonds), and 0.050 (triangles). Top panel is the result at m s = 0.08 and 
the bottom is in m s = 0.10. 



C. Non-conserved current 

The discretization effect may also come from the non-conserved local current J x ° c in (fT3]) , 
which is represented by the term A^ U (Q) in (11 Op . 

We estimate its leading contribution A^ U (Q) = (S l _ lu Q 2 —Q^ t Q u )Q 2 ,F^ 1 (Q) in the expansion 
f TT"3]) in terms of small aQ^, by solving linear equations ( ITU]) for different sets of [i and v. 
We find that the maximum magnitude of (S^Q 2 — Q f _ l Q u )Ql^ii +A {Q) is less than 1% of 
U V+A {Q) in the fit range 0.463 < {aQ) 2 < 0.994. 

The contribution of A^(Q) may also be estimated by looking at a difference between 
VPFs obtained with ji = v (diag) and with /i ^ v (offd) components, i.e. 

nf g (Q) = (J?J l ° c )(Q)/(Q 2 - Q»Q»), (15) 
n°/ d (Q) = (J?J l ° c )(Q)/(-QM, (16) 

respectively. Figure [6] shows Uy as (Q) — ily fd (Q) as a function of (aQ) 2 . The maximum 
magnitude of the difference in the range 0.463 < (aQ) 2 < 0.994 is 0.003, which is the same 
order as the estimate from F^ +A (Q). 

Adding or subtracting this amount of systematic effect from the lattice data, we repeat 
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FIG. 7: C (Q 2 , fx 2 , a s ) — c + l/(2ir 2 )ln(aQ) 2 as a function of (aQ) 2 . The perturbative results at 
two-loop (solid), three- loop (dashed), and four- loop (dotted) calculations are shown. The logarithm 
at the leading order is subtracted in order to enhance their small differences. 



the whole analysis to estimate the systematic error due to the Lorentz (or WT) violating 
terms, which gives ±0.0002 in a^(Mz)- 

D. Quark condensate and Z m 

The uncertainty due to the quark condensate is estimated as ±0.0001 for ai (Mz) by 
varying the input value from —[0.220 GeV] 3 to —[0.265 GeV] 3 , which correspond to the lower 



and upper limits of the estimate of (qq) in 23]. 



The uncertainty due to Z m is also estimated as ±0.0001 for a^(Mz) by varying Z m 
within its estimated error (from 0.777 to 0.832). 

E. Perturbative expansion 

The truncation effect of the perturbative expansion can be estimated by comparing the 
results with different orders of the perturbative expansion. Fortunately, the four-loop cal- 



culation is known for Co(Q 2 , fi 2 , a s ) [28|, [29[, and we can explicitly estimate the effect of 
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A comparison of two-, three- and four-loop calculations of Cq(Q 2 , fi 2 ,a s ) is shown in 
Figure [3 They correspond to 0(a s ), 0(a 2 ), and 0(a%) calculations, respectively. We 
observe that the difference between three-loop and four-loop is of order of 0.0001 for Uj(Q), 
which is much smaller than other systematic effect. 

Strictly speaking, the smallness of the four-loop contribution does not guarantee that 
the unknown higher orders are even smaller. We therefore attempted to fit the data with a 
formula including unknown O(aj) term c^aj(Q). The fit gives ~ O(10) with a shift of 
resulting af\Mz) by +0.0003. We therefore put a conservative systematic error from the 
truncation of the perturbation series as ±0.0003. This can only be reduced by including the 
data at higher Q 2 values, which needs finer lattice spacings. 

F. 1/Q 2 expansion 

As previously discussed, the size of neglected 1/Q 6 terms in the OPE formula is at most 
0.001 for Uj(Q) at the lower end of the fitting range (aQ)^ in = 0.463 (see Figure H]). This 
is less than 1/5 of the estimated discretization effect discussed above. We therefore expect 
that the impact on af'(M^) is smaller than 0.0001. 

G. Charm and bottom quark mass 

The uncertainty of charm and bottom quark masses, m c and nit, used in a perturbative 
matching of ai 3 ^ onto ai 5 ^ is +0.0001 and —0.0003, which are the maximum and mini- 
mum values when mb tC are changed within la in the analysis. The input values m c (m c ) = 
1.27(i°[) GeV and m b (rh b ) = 4.20(+^) GeV are taken from Q. 

H. Lattice spacing 

The uncertainty of the lattice spacing is the largest source of the systematic error. 

Our main result is quoted with the Sommer scale tq = 0.49 fm as an input, with which 
we obtain a~ l = 1.83(1) GeV. This quantity is convenient because its numerical calculation 
is very precise and also because its sea quark mass dependence is mild. On the other hand, 
r does not have a direct relation to any physical observables and one has to resort to some 
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model to fix the central value. For this reason, one prefers other physical quantities to set 
the scale. 

One possible candidate is the pion decay constant /„., which is also precisely measurable 
on the lattice at unphysical values of sea quark masses. The problem for this quantity 
is that it may have rather non-trivial sea quark mass dependence as predicted by Chiral 



Perturbation Theory (ChP' 
obtain a" 1 = 1.97(4) GeV 



Y Using the next-to-next-to-leading order formula in ChPT, we 



30j ] . In this analysis we observed a rather large dependence on 



the sea quark mass and more importantly a strong curvature that bends the extrapolation 
to the lower value of af n , thus the higher value of a -1 . We therefore need more careful 
analysis on possible systematic errors in this determination. 

One of the most popular quantities to set the scale in recent lattice calculations is the Q 
baryon mass. Since the Q baryon is made of three strange quarks, the dependence on up 
and down quark masses only comes from quark loop effect, which is expected to be small. 
A possible problem is that its determination has to be combined with the strange quark 
mass determination, which is non-trivial. In addition, the finite volume effect could be more 
important for baryons. Our result is a' - 1 = 1.76(8) (tl) GeV, with the second error being 
our estimate of the finite volume effect, that is set by calculations on a larger volume lattice 
(24 3 x 48) but at limited values of sea quark mass. 

Since each determination has its own advantage and disadvantage, we decided to take 
ro as our central value and others (/^ and Mq) to estimate the systematic uncertainties. 
The shift of a s (Mz) due to the choice of f n and Mq is +0.0013 and —0.0010, respectively, 
which we quote as the systematic error from the scale setting. 

This uncertainty also affects the matching points m c and mj, which is included in the 
above error band. 



I. Final result 

Table |T] shows a summary of our estimate of the systematic errors in our determination 
of a s (Mz)- We quote final result of the strong coupling constant at the Z boson mass scale 
as 

af\M z )= 0.1181(3)0. (17) 
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Sources 


Estimated error in a^{M z ) 


Uncorrelated fit 


±0.0003 


Lattice artifact {0(a?) effect) 


+0.0003 




±0.0002 


Quark condensate 


±0.0001 




±0.0001 


Perturbative expansion 


±0.0003 


1/Q 2 expansion 


< 0.0001 


m c>b 


+0.0001 


-0.0003 


Lattice spacing 


+0.0013 
-0.0010 


Total (in quadrature) 


+0.0014 
-0.0012 



TABLE I: Summary of systematic errors in a s (Mz). 



Here, the first error is statistical error and the second is a sum of the various systematic 
errors in quadrature. 

This result is consistent with other recent lattice QCD results 0.1174(12) {q], 0.1183(8) 



various high-energy experiments |l|, |2[ (updated on-line in 2010) 
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3l|, 0.1192(11) [32J, 0.1205(8) (5) (l^)[33j, and with the world average 0.1184(7) including 



V. CONCLUSIONS 



Determination of the strong coupling constant a s (fi) may be achieved through a pertur- 
bative expansion of any physical quantity in terms of a s (fi) at a given scale fi. Experimental 
determination typically uses a perturbative amplitude of quarks at high energy. Comparison 
with the lattice QCD calculation provides a highly non-trivial test of QCD, as lattice uses 
low-energy hadron spectrum or matrix elements to set the scale. 

In the lattice calculation, there are variety of choices for the quantity to be expanded 
in a s (fi). In order to achieve good enough accuracy, the perturbative expansion must be 
known to higher orders or evaluated at very high energies. The latter may be achieved 
by calculating the scaling towards the high energy regime non-perturbatively using the so- 
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called step-scaling technique (see, e.g. [5|, |33|). The former is numerically less intensive but 
requires analytic perturbative calculations beyond one-loop level. 

This work demonstrates that the vacuum polarization function can be used for the precise 
determination of a s . The important points are (i) the perturbative expansion can be done in 
the continuum theory and is known to O(a^), and (ii) the non-perturbative lattice calculation 
with controlled systematic errors is possible. The discretization error was a concern as the 
large Q 2 points have to be calculated, but it turned out to be under control with currently 
available lattice setups by careful estimates of systematic effects. 

The use of the overlap fermion is certainly desirable as the massless limit of quarks is 
uniquely defined and the use of the continuum OPE is justified. With the lattice fermions 
that violates chiral symmetry, one expects dangerous terms such as ma _3 /Q 4 , whose nu- 
merical impact has to be carefully studied. 

Extension of this work is straightforward. Since the largest uncertainty comes from 
the scale determination, a consistent determination of the lattice scale with various low- 
energy inputs is necessary in order to significantly improve the accuracy. This requires 
extensive simulations at larger volumes, smaller quark masses and smaller lattice spacings. 
The discretization effect in VPF considered in this work will also be significantly reduced 
by going to finer lattice that will become available within a few years. 
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